Plots

rm(list=ls());gc()

##          used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 411635 22.0     851872 45.5   648454 34.7
## Vcells 756937  5.8    8388608 64.0  1699667 13.0

load("Definitive_models_2021.RData")
#xaringan::inf_mr()

#arriba puse algunas opciones para que por defecto escondiera el código
#también cargue algunos estilo .css para que el texto me apareciera justificado, entre otras cosas.
local({r <- getOption("repos")
       r["CRAN"] <- "http://cran.r-project.org" 
       options(repos=r)
})

`%>%` <- magrittr::`%>%`
copy_names <- function(x,row.names=FALSE,col.names=TRUE,dec=",",...) {
  if(class(ungroup(x))[1]=="tbl_df"){
        if(options()$OutDec=="."){
            options(OutDec = dec)
            write.table(format(data.frame(x)),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
            options(OutDec = ".")
          return(x)
        } else {
            options(OutDec = ",")
            write.table(format(data.frame(x)),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
            options(OutDec = ",")
          return(x)    
        }
  } else {
        if(options()$OutDec=="."){
            options(OutDec = dec)
            write.table(format(x),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
            options(OutDec = ".")
          return(x)
        } else {
            options(OutDec = ",")
            write.table(format(x),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
            options(OutDec = ",")
          return(x)       
  }
 }
}  

unlink('Plots 2021-01-18_cache', recursive = TRUE)
if(!require(pacman)){install.packages("pacman")}

pacman::p_unlock(lib.loc = .libPaths()) #para no tener problemas reinstalando paquetes
knitr::opts_chunk$set(
    echo = TRUE,
    message = FALSE,
    warning = FALSE
)
#dejo los paquetes estadísticos que voy a utilizar

if(!require(plotly)){install.packages("plotly")}
if(!require(lubridate)){install.packages("lubridate")}
if(!require(htmlwidgets)){install.packages("htmlwidgets")}
if(!require(tidyverse)){install.packages("tidyverse")}
if(!require(gganimate)){install.packages("gganimate")}
if(!require(readr)){install.packages("readr")}
if(!require(stringr)){install.packages("stringr")}
if(!require(data.table)){install.packages("data.table")}
if(!require(DT)){install.packages("DT")}
if(!require(ggplot2)){install.packages("ggplot2")}
if(!require(lattice)){install.packages("lattice")}
if(!require(forecast)){install.packages("forecast")}
if(!require(zoo)){install.packages("zoo")}
if(!require(panelView)){install.packages("panelView")}
if(!require(janitor)){install.packages("janitor")}
if(!require(rjson)){install.packages("rjson")}
if(!require(estimatr)){install.packages("estimatr")} 
if(!require(CausalImpact)){install.packages("CausalImpact")}
if(!require(textreg)){install.packages("textreg")}
if(!require(sjPlot)){install.packages("sjPlot")}
if(!require(foreign)){install.packages("foreign")}
if(!require(tsModel)){install.packages("tsModel")}
if(!require(lmtest)){install.packages("lmtest")}
if(!require(Epi)){install.packages("Epi")}
if(!require(splines)){install.packages("splines")}
if(!require(vcd)){install.packages("vcd")}
if(!require(astsa)){install.packages("astsa")}
if(!require(forecast)){install.packages("forecast")}
if(!require(MASS)){install.packages("MASS")}
if(!require(ggsci)){install.packages("ggsci")}
if(!require(Hmisc)){install.packages("Hmisc")}
if(!require(compareGroups)){install.packages("compareGroups")}
if(!require(dplyr)){install.packages("dplyr")}
if(!require(ggforce)){install.packages("ggforce")}
if(!require(imputeTS)){install.packages("imputeTS")}
if(!require(doParallel)){install.packages("doParallel")}
if(!require(SCtools)){install.packages("SCtools")}
if(!require(MSCMT)){install.packages("MSCMT")}
# Calculate the number of cores
no_cores <- detectCores() - 1
cl<-makeCluster(no_cores)
registerDoParallel(cl)

Sys.setlocale(category = "LC_ALL", locale = "english")

## [1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"

0.1 Final Plots

library(cowplot)

#horizontal
height_y<-16
height_x<-16
size_title<-18
line_size<-1.2 #.8


Sys.setlocale(category = "LC_ALL", locale = "english")

## [1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"

hosp_trauma_plot_red<-
plot(impact3d_hosp_trauma, "original")$data %>% 
    ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=format(as.Date(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])),"%Y %b"))+
  ylim(c(0,150))+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("c) Trauma Hospitalizations")

hosp_resp_plot_red<-
plot(impact3d1_hosp_resp, "original")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=format(as.Date(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])),"%Y %b"))+
  ylim(c(0,150))+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("d) Respiratory Hospitalizations")

cons_trauma_plot_red<-
plot(impact3d1_cons_trauma, "original")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=format(as.Date(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])),"%Y %b"))+
  ylim(c(0,1600))+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("a) Trauma Consultations")

cons_resp_plot_red<-
plot(impact3d_cons_resp, "original")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=format(as.Date(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])),"%Y %b"))+
  ylim(c(0,1600))+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("b) Respiratory Consultations")

rate_trauma_plot_red<-
plot(impact3d_cons_resp, "original")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=format(as.Date(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])),"%Y %b"))+
  ylim(c(0,400))+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("e) Trauma Hospitalizations per 1,000 consultations")

rate_resp_plot_red<-
plot(impact3d_ratio_resp, "original")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="red") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=format(as.Date(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])),"%Y %b"))+
  ylim(c(0,400))+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("f) Respiratory Hospitalizations per 1,000 consultations")


#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

hosp_trauma_plot<-
plot(impact3d_hosp_trauma, "original")+
    xlab("")+
    #ylab("Hospitalizations")+
    scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=as.character(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])))+
    ylim(c(0,150))+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("c) Trauma Hospitalizations")
hosp_resp_plot<-
plot(impact3d1_hosp_resp, "original")+
    xlab("")+
    #ylab("Hospitalizations")+
    scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=as.character(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])))+
    ylim(c(0,150))+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("d) Respiratory Hospitalizations")
cons_trauma_plot<-
plot(impact3d1_cons_trauma, "original")+
    xlab("")+
    #ylab("Consultations")+
    scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=as.character(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])))+
    ylim(c(0,1600))+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("a) Trauma Consultations")
cons_resp_plot<-
plot(impact3d_cons_resp, "original")+
    xlab("")+
    #ylab("Consultations")+
    scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=as.character(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])))+
    ylim(c(0,1600))+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("b) Respiratory Consultations")
rate_trauma_plot<-
plot(impact3d_ratio, "original")+
    xlab("")+
    #ylab("Rate of Hospitalizations per Consultations")+
    ylim(c(0,400))+
    scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=as.character(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])))+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("e) Trauma Hospitalizations per 1,000 consultations")
rate_resp_plot<-
plot(impact3d_ratio_resp, "original")+
    xlab("")+
    #ylab("Rate of Hospitalizations per Consultations")+
    ylim(c(0,400))+
    scale_x_continuous(breaks=(c(seq(1,nrow(data15a64_rn),24),262)),
                       labels=as.character(unlist(data15a64_rn[c(seq(1,nrow(data15a64_rn),24),262),"date"])))+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("f) Respiratory Hospitalizations per 1,000 consultations")

plot1<-
plot_grid(
  cons_trauma_plot+theme(axis.text.x = element_blank()),
  cons_resp_plot+theme(axis.text.x = element_blank()),
  hosp_trauma_plot+theme(axis.text.x = element_blank()),
  hosp_resp_plot+theme(axis.text.x = element_blank()),
  rate_trauma_plot+theme(axis.text.x = element_blank()),
  rate_resp_plot,
  #get_legend(t14_trend_leg + theme(legend.position="right",
  #                                 legend.text = element_text(size = 15))),
  nrow = 6,  rel_heights = c(rep(1,5),1.52), align = "v"
#, rel_widths = c(1, 1)
)

#vertical
height_y<-13
height_x<-13
size_title<-14

plot1b<-
plot_grid(
  cons_trauma_plot+theme(axis.text.x = element_blank()),
  cons_resp_plot+theme(axis.text.x = element_blank()),
  hosp_trauma_plot+theme(axis.text.x = element_blank()),
  hosp_resp_plot+theme(axis.text.x = element_blank()),
  rate_trauma_plot,
  rate_resp_plot,#,
  nrow = 3, rel_widths = c(1, 1), rel_heights = c(rep(1,3),1.7)
)

plot1b_red<-
plot_grid(
  cons_trauma_plot_red+theme(axis.text.x = element_blank()),
  cons_resp_plot_red+theme(axis.text.x = element_blank()),
  hosp_trauma_plot_red+theme(axis.text.x = element_blank()),
  hosp_resp_plot_red+theme(axis.text.x = element_blank()),
  rate_trauma_plot_red,
  rate_resp_plot_red,#,
  nrow = 3, rel_widths = c(1, 1), rel_heights = c(rep(1,2),1.25)
)


ggdraw(plot1b_red)+
   theme(plot.background = element_rect(fill=NA, color = NA))

Figure 11. Estimated Trends of the Outcomes

dont_show=1
if(dont_show==0){
  pdf("_Fig1.pdf", height = 14, width = 23)
  ggdraw(plot1)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  jpeg("_Fig1.jpg", height = 14, width = 23, units = 'in', res = 600)
  ggdraw(plot1)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  pdf("_Fig1b.pdf", height = 19.5, width = 15.3)
  ggdraw(plot1b)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  jpeg("_Fig1b.jpg", height = 19.5, width = 15.3, units = 'in', res = 600)
  ggdraw(plot1b)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}

if(dont_show==0){
  pdf("_Fig1c.pdf", height = 14, width = 23)
  ggdraw(plot1b_red)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  jpeg("_Fig1c.jpg", height = 14, width = 23, units = 'in', res = 600)
  ggdraw(plot1b_red)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
#add_sub(, "Note. The first panel shows the data and a counterfactual prediction for the post-treatment period (Blue dashed line);\nBlue area= Prediction intervals.", vpadding=grid::unit(0, "lines"), y = .55, x = 0.03, hjust = 0,size=9)
#plot2 <- cowplot::ggdraw(grid.arrange(p14, p21,p212,p32,p42,p52,p57, ncol = 2, nrow = 4)) + 
  # same plot.background should be in the theme of p1 and p2 as mentioned above
#  theme(plot.background = element_rect(fill=NA, color = NA))

#horizontal
height_y<-16
height_x<-16
size_title<-18
line_size<-.8

hosp_trauma_plot2_red<-
plot(impact3d_hosp_trauma, "pointwise")$data %>% 
    ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("c) Trauma Hospitalizations")

hosp_resp_plot2_red<-
plot(impact3d1_hosp_resp, "pointwise")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("d) Respiratory Hospitalizations")

cons_trauma_plot2_red<-
plot(impact3d1_cons_trauma, "pointwise")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("a) Trauma Consultations")

cons_resp_plot2_red<-
plot(impact3d_cons_resp, "pointwise")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("b) Respiratory Consultations")

rate_trauma_plot2_red<-
plot(impact3d_ratio, "pointwise")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("e) Trauma Hospitalizations per 1,000 consultations")

rate_resp_plot2_red<-
plot(impact3d_ratio_resp, "pointwise")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("f) Respiratory Hospitalizations per 1,000 consultations")

vector_breaks_plot2b_red<- seq(238,nrow(data15a64_rn),3)
Sys.setlocale(category = "LC_ALL", locale = "english")

## [1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"

vector_dates_plot2b_red<-as.character(format(as.Date(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),"%b %d"))

plot2b_red<-
plot_grid(
  cons_trauma_plot2_red+scale_x_continuous(breaks=vector_breaks_plot2b_red,labels=vector_dates_plot2b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-800,300,200), limits = c(-800,300)),
  cons_resp_plot2_red+scale_x_continuous(breaks=vector_breaks_plot2b_red,labels=vector_dates_plot2b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-800,300,200), limits = c(-800,300)),
  hosp_trauma_plot2_red+scale_x_continuous(breaks=vector_breaks_plot2b_red,labels=vector_dates_plot2b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-75,75,25), limits = c(-75,75)),
  hosp_resp_plot2_red+scale_x_continuous(breaks=vector_breaks_plot2b_red,labels=vector_dates_plot2b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-75,75,25), limits = c(-75,75)),
  rate_trauma_plot2_red+scale_x_continuous(breaks=vector_breaks_plot2b_red,labels=vector_dates_plot2b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-400,400,200), limits = c(-400,400)),
  rate_resp_plot2_red+scale_x_continuous(breaks=vector_breaks_plot2b_red,labels=vector_dates_plot2b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-400,400,200), limits = c(-400,400)),
  #get_legend(t14_trend_leg + theme(legend.position="right",
  #                                 legend.text = element_text(size = 15))),
    nrow = 3, rel_widths = c(.95, .95), rel_heights = c(rep(1,2),1.3), scale = c(.95, .95, .95,.95,.95,.95)
)+
  draw_label("Weekly difference change in consultations/hospitalizations",x=0.005, y=.5, angle=90, size = 17)

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

#vertical
height_y<-13
height_x<-13
size_title<-14

hosp_trauma_plot2<-
plot(impact3d_hosp_trauma, "pointwise")+
    xlab("")+
    #ylab("Hospitalizations")+
    xlim(238,262)+
    theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("c) Trauma Hospitalizations")+
  theme(plot.title = element_text(size=size_title))
hosp_resp_plot2<-
plot(impact3d1_hosp_resp, "pointwise")+
    xlab("")+
    #ylab("Hospitalizations")+
    xlim(238,262)+
    theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("d) Respiratory Hospitalizations")+
  theme(plot.title = element_text(size=size_title))
cons_trauma_plot2<-
plot(impact3d1_cons_trauma, "pointwise")+
    xlab("")+
    #ylab("Consultations")+
    xlim(238,262)+
    theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("a) Trauma Consultations")+
  theme(plot.title = element_text(size=size_title))
cons_resp_plot2<-
plot(impact3d_cons_resp, "pointwise")+
    xlab("")+
    #ylab("Consultations")+
    xlim(238,262)+
    theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=12))+
  ggtitle("b) Respiratory Consultations")+
  theme(plot.title = element_text(size=size_title))
rate_trauma_plot2<-
plot(impact3d_ratio, "pointwise")+
    xlab("")+
    #ylab("Rate of Hospitalizations per Consultations")+
    xlim(238,262)+
    theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("e) Trauma Hospitalizations per 1,000 consultations")+
  theme(plot.title = element_text(size=size_title))
rate_resp_plot2<-
plot(impact3d_ratio_resp, "pointwise")+
    xlab("")+
    #ylab("Rate of Hospitalizations per Consultations")+
    xlim(238,262)+
    theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("f) Respiratory Hospitalizations per 1,000 consultations")+
  theme(plot.title = element_text(size=size_title))

Sys.setlocale(category = "LC_ALL", locale = "english")

## [1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"

plot2<-
plot_grid(
  cons_trauma_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=format(as.Date(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),"%b %d"),limits=c(238,262))+scale_y_continuous(breaks=seq(-800,300,200), limits = c(-800,300)),
  cons_resp_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=format(as.Date(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),"%b %d"),limits=c(238,262))+scale_y_continuous(breaks=seq(-800,300,200), limits = c(-800,300)),
  hosp_trauma_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=format(as.Date(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),"%b %d"),limits=c(238,262))+scale_y_continuous(breaks=seq(-75,75,25), limits = c(-75,75)),
  hosp_resp_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=format(as.Date(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),"%b %d"),limits=c(238,262))+scale_y_continuous(breaks=seq(-75,75,25), limits = c(-75,75)),
  rate_trauma_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=format(as.Date(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),"%b %d"),limits=c(238,262))+scale_y_continuous(breaks=seq(-100,400,100), limits = c(-100,400)),
  rate_resp_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=format(as.Date(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),"%b %d"),limits=c(238,262))+scale_y_continuous(breaks=seq(-100,400,100), limits = c(-100,400)),
  #get_legend(t14_trend_leg + theme(legend.position="right",
  #                                 legend.text = element_text(size = 15))),
  nrow = 3, rel_widths = c(1, 1)
)+
  draw_label("Weekly difference change in consultations/hospitalizations",x=0.01, y=.5, angle=90, size = 14)

#horizontal
height_y<-15
height_x<-15
size_title<-16


Sys.setlocale(category = "LC_ALL", locale = "english")

## [1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"

plot2b<-
plot_grid(
  cons_trauma_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=format(as.Date(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),"%b %d"),limits=c(238,262))+scale_y_continuous(breaks=seq(-800,300,200), limits = c(-800,300)),
  cons_resp_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=format(as.Date(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),"%b %d"),limits=c(238,262))+scale_y_continuous(breaks=seq(-800,300,200), limits = c(-800,300)),
  hosp_trauma_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=format(as.Date(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),"%b %d"),limits=c(238,262))+scale_y_continuous(breaks=seq(-75,75,25), limits = c(-75,75)),
  hosp_resp_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=format(as.Date(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),"%b %d"),limits=c(238,262))+scale_y_continuous(breaks=seq(-75,75,25), limits = c(-75,75)),
  rate_trauma_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=format(as.Date(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),"%b %d"),limits=c(238,262))+scale_y_continuous(breaks=seq(-100,400,100), limits = c(-100,400)),
  rate_resp_plot2+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=format(as.Date(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),"%b %d"),limits=c(238,262))+scale_y_continuous(breaks=seq(-100,400,100), limits = c(-100,400)),
  #get_legend(t14_trend_leg + theme(legend.position="right",
  #                                 legend.text = element_text(size = 15))),
  nrow = 3, rel_widths = c(1, 1)
)+
  draw_label("Weekly difference change in consultations/hospitalizations",x=0.01, y=.5, angle=90, size = 14)
 #draw_figure_label(label = "\n\nFigure 1", angle = -90, colour = "black", position="bottom.left", size=20)
  
    ggdraw(plot2b_red)+
     theme(plot.background = element_rect(fill=NA, color = NA))

Figure 12. Estimated Trends of the Outcomes

#add_sub(, "Note. Difference between observed data and counterfactual predictions;\nBlue area= Prediction intervals.", 
               #vpadding=grid::unit(0, "lines"),
                #   y = .55, x = 0.03, hjust = 0,size=9)

if(dont_show==0){
  pdf("___Fig2.pdf", width = 15.3, height = 19.5)
  ggdraw(plot2)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  jpeg("_Fig2.jpg", width = 15.3, height = 19.5, units = 'in', res = 600)
  ggdraw(plot2)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  pdf("_Fig2b.pdf", height = 14, width = 23)
  ggdraw(plot2b)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  jpeg("_Fig2b.jpg", height = 14, width = 23, units = 'in', res = 600)
  ggdraw(plot2b)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  pdf("_Fig2c.pdf", height = 14, width = 23)
  ggdraw(plot2b_red)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  jpeg("_Fig2c.jpg", height = 14, width = 23, units = 'in', res = 600)
  ggdraw(plot2b_red)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}

#horizontal
height_y<-16
height_x<-16
size_title<-18
line_size<-.8

hosp_trauma_plot3_red<-
plot(impact3d_hosp_trauma, "cumulative")$data %>% 
    ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("c) Trauma Hospitalizations")

hosp_resp_plot3_red<-
plot(impact3d1_hosp_resp, "cumulative")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("d) Respiratory Hospitalizations")

cons_trauma_plot3_red<-
plot(impact3d1_cons_trauma, "cumulative")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("a) Trauma Consultations")

cons_resp_plot3_red<-
plot(impact3d_cons_resp, "cumulative")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_blank(),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("b) Respiratory Consultations")

rate_trauma_plot3_red<-
plot(impact3d_ratio, "cumulative")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("e) Trauma Hospitalizations per 1,000 consultations")

rate_resp_plot3_red<-
plot(impact3d_ratio_resp, "cumulative")$data %>% 
      ggplot(aes(x=time))+
      geom_line(aes(y=mean), size=line_size, linetype="longdash", color="darkblue") +
      geom_line(aes(y=response), size=line_size, color="black")+
  geom_ribbon(aes(ymin=lower,ymax=upper),fill="steelblue", alpha = 0.35)+
  theme_sjplot()+
  xlim(238,262)+
  geom_hline(yintercept = 0, col = "black", lty = 3, size=1)+
  geom_vline(xintercept = length(data15a64_rn$did[data15a64_rn$did==0]), col = "gray50", lty = 2, size=1)+
      theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=.5,size=height_x),
     # theme(axis.text.x = element_blank(),#element_text(angle = 90, vjust = 0.5, hjust=.5,size=9),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_blank(),
          axis.title.x = element_blank(),
          plot.title = element_text(size=size_title),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
    ggtitle("f) Respiratory Hospitalizations per 1,000 consultations")

vector_breaks_plot3b_red<- seq(238,nrow(data15a64_rn),3)
Sys.setlocale(category = "LC_ALL", locale = "english")

## [1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"

vector_dates_plot3b_red<-as.character(format(as.Date(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),"%b %d"))

plot3b_red<-
plot_grid(
  cons_trauma_plot3_red+scale_x_continuous(breaks=vector_breaks_plot3b_red,labels=vector_dates_plot3b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-5000,2000,1000), limits = c(-5000,2000)),
  cons_resp_plot3_red+scale_x_continuous(breaks=vector_breaks_plot3b_red,labels=vector_dates_plot3b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-5000,2000,1000), limits = c(-5000,2000)),
  hosp_trauma_plot3_red+scale_x_continuous(breaks=vector_breaks_plot3b_red,labels=vector_dates_plot3b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-250,250,100), limits = c(-250,250)),
  hosp_resp_plot3_red+scale_x_continuous(breaks=vector_breaks_plot3b_red,labels=vector_dates_plot3b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-250,250,100), limits = c(-250,250)),
  rate_trauma_plot3_red+scale_x_continuous(breaks=vector_breaks_plot3b_red,labels=vector_dates_plot3b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-100,1500,250), limits = c(-100,1500)),
  rate_resp_plot3_red+scale_x_continuous(breaks=vector_breaks_plot3b_red,labels=vector_dates_plot3b_red,limits=c(238,262))+scale_y_continuous(breaks=seq(-100,1500,250), limits = c(-100,1500)),
  #get_legend(t14_trend_leg + theme(legend.position="right",
  #                                 legend.text = element_text(size = 15))),
    nrow = 3, rel_widths = c(1, 1), rel_heights = c(rep(1,2),1.25), scale = c(.95, .95, .95,.95,.95,.95))+
 draw_label("Weekly difference change in consultations/hospitalizations",x=0.005, y=.5, angle=90, size = 17)

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

#
#vertical
height_y<-13
height_x<-13
size_title<-14

#horizontal

hosp_trauma_plot3<-
plot(impact3d_hosp_trauma, "cumulative")+
    xlab("")+
   # ylab("Hospitalizations")+
      xlim(238,262)+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("c) Trauma Hospitalizations")+
  theme(plot.title = element_text(size=size_title))
hosp_resp_plot3<-
plot(impact3d1_hosp_resp, "cumulative")+
    xlab("")+
   # ylab("Hospitalizations")+
      xlim(238,262)+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("d) Respiratory Hospitalizations")+
  theme(plot.title = element_text(size=size_title))
cons_trauma_plot3<-
plot(impact3d1_cons_trauma, "cumulative")+
    xlab("")+
   # ylab("Consultations")+
    xlim(238,262)+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("a) Trauma Consultations")+
  theme(plot.title = element_text(size=size_title))
cons_resp_plot3<-
plot(impact3d_cons_resp, "cumulative")+
    xlab("")+
   # ylab("Consultations")+
    xlim(238,262)+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("b) Respiratory Consultations")+
  theme(plot.title = element_text(size=size_title))
rate_trauma_plot3<-
plot(impact3d_ratio, "cumulative")+
    xlab("")+
   # ylab("Rate of Hospitalizations per Consultations")+
    xlim(238,262)+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("e) Trauma Hospitalizations per 1,000 consultations")+
  theme(plot.title = element_text(size=size_title))
rate_resp_plot3<-
plot(impact3d_ratio_resp, "cumulative")+
    xlab("")+
   # ylab("Rate of Hospitalizations per Consultations")+
    xlim(238,262)+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=.5,size=height_x),
          axis.text.y = element_text(size=height_y),
          axis.title.y = element_text(size=9),
          axis.title.x = element_blank(),
          plot.caption = element_text(hjust = 0, face= "italic",size=9))+
  ggtitle("f) Respiratory Hospitalizations per 1,000 consultations")+
  theme(plot.title = element_text(size=size_title))

plot3<-
plot_grid(
  cons_trauma_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-4000,2000,1000), limits = c(-4000,2000)),
  cons_resp_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-4000,2000,1000), limits = c(-4000,2000)),
  hosp_trauma_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-250,250,100), limits = c(-250,250)), 
  hosp_resp_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-250,250,100), limits = c(-250,250)),
  rate_trauma_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-100,1500,250), limits = c(-100,1500)),
  rate_resp_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-100,1500,250), limits = c(-100,1500)),#,
  #get_legend(t14_trend_leg + theme(legend.position="right",
  #                                 legend.text = element_text(size = 15))),
  nrow = 3, rel_widths = c(1, 1)
)

height_y<-15
height_x<-15
size_title<-16

plot3b<-
plot_grid(
  cons_trauma_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-4000,2000,1000), limits = c(-4000,2000)),
  cons_resp_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-4000,2000,1000), limits = c(-4000,2000)),
  hosp_trauma_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-250,250,100), limits = c(-250,250)), 
  hosp_resp_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-250,250,100), limits = c(-250,250)),
  rate_trauma_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-100,1500,250), limits = c(-100,1500)),
  rate_resp_plot3+scale_x_continuous(breaks=seq(238,nrow(data15a64_rn),3),labels=as.character(unlist(data15a64_rn[c(seq(238,nrow(data15a64_rn),3)),"date"])),limits=c(238,262))+scale_y_continuous(breaks=seq(-100,1500,250), limits = c(-100,1500)),#,
  #get_legend(t14_trend_leg + theme(legend.position="right",
  #                                 legend.text = element_text(size = 15))),
  nrow = 3, rel_widths = c(1, 1)
)
#add_sub(, "Note. Added pointwise contributions from the second panel;\nBlue area= Cumulative intervals;\nBlue area= Prediction intervals.", 
               #vpadding=grid::unit(0, "lines"),
                #   y = .55, x = 0.03, hjust = 0,size=9))+
ggdraw(plot3b_red)+
   theme(plot.background = element_rect(fill=NA, color = NA))

Figure 13. Estimated Trends of the Outcomes

if(dont_show==0){
  pdf("___Fig3.pdf", width = 15.3, height = 19.5)
  ggdraw(plot3)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  jpeg("_Fig3.jpg", width = 15.3, height = 19.5, units = 'in', res = 600)
  ggdraw(plot3)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  pdf("_Fig3b.pdf", height = 14, width = 23)
  ggdraw(plot3b)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  jpeg("_Fig3b.jpg", height = 14, width = 23, units = 'in', res = 600)
  ggdraw(plot3b)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  pdf("_Fig3c.pdf", height = 14, width = 23)
  ggdraw(plot3b_red)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}
if(dont_show==0){
  jpeg("_Fig3c.jpg", height = 14, width = 23, units = 'in', res = 600)
  ggdraw(plot3b_red)+
     theme(plot.background = element_rect(fill=NA, color = NA))
  dev.off()
}